rm(list=ls(all=TRUE)) 
results <- read.csv("Data/Urchin_ClimateChange_Data.csv", header=T, na="NA")

library("lsmeans") #post-hoc tests
library("effects") #plot effects of modeling
library("lmtest") #linear mixed modeling
library("lme4") #linear mixed modeling
library("lmerTest") #calculate p-values in models
library("emmeans") #post-hoc tests
library("cowplot") #arrange plots
library("MuMIn")
library("car") #levenes tests, Anova, qqPlots, 
library("tidyverse")
library("stats")
library("plyr")  #splitting, applying, and combining data
library("dplyr")
library("onewaytests") #allows for Welch test with unequal variances (spine lengths)
library("stats")
library("effsize")
library("FSA")
library("rcompanion")
library("broom") #turns model output into data frames
library("plotrix") #Calculates standard error of the mean

Analysis approach:

  1. Build and run model
  2. Check for normality of residuals
  3. Check for homogeniety of variance of residuals
  4. Look at model summary
  5. Run anova to check for significance of main effects
  6. Run post-hoc to test hypotheses

Data Assembly and Results:

Diameter Means

Day -24

Assemble data for diameters of initial collection of urchins from hatchery (see code)
##Initial Collection = Day -24
InitialDiameter <-
  #seperate out the initial sizes of urchins on day -24, when desired conditions were met.
  results %>% 
    select("Day", "Diam1", "Diam2", "TankID", "Temperature", "pH", "Treatment") %>%
    filter(TankID != "8t" ) %>% 
      #remove urchin number 8 which died halfway through
    drop_na(Diam1, Diam2) %>% 
    gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>% 
      #gather Diam1 and Diam2 into single column in order to make number of levels of each grouping factor < number of observations (as per error message previously recieved)
    filter(Day == "-24") %>% 
      #sizes on day -24  
    select("TankID", "Diameter", "Temperature", "pH", "Treatment")


InitialMean <- mean(InitialDiameter$Diameter)
InitialSE <- std.error(InitialDiameter$Diameter)
Results of linear mixed model effect on diameter (day -24)
InitialMod <- lmer(Diameter ~ Temperature*pH + (1|TankID), data=InitialDiameter)
anova(InitialMod, type = 2)
## Type II Analysis of Variance Table with Satterthwaite's method
##                   Sum Sq   Mean Sq NumDF DenDF F value Pr(>F)
## Temperature    0.0005035 0.0005035     1    19  0.1144 0.7389
## pH             0.0102269 0.0102269     1    19  2.3232 0.1439
## Temperature:pH 0.0003840 0.0003840     1    19  0.0872 0.7709
Check assumptions of diameter model (day -24)
#check assumptions
# 1. Normality of residuals
qqPlot(residuals(InitialMod)) 

## [1] 27  4
shapiro.test(residuals(InitialMod)) #no pass
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(InitialMod)
## W = 0.52108, p-value = 4.878e-11
hist(residuals(InitialMod)) #eh

# 2. Equal variances
leveneTest(residuals(InitialMod)~InitialDiameter$Treatment) #No pass
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  3  2.9592 0.04307 *
##       42                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(fitted(InitialMod),resid(InitialMod,type="pearson"),col="blue")

Day -13

Assemble Data for diameters of acclimated urchins before conditioning (see code)
##After acclimating before ramp up = Day-13

AcclimatedDiameter <-
  #seperate out the initial sizes of urchins on day -24, when desired conditions were met.
  results %>% 
    select("Day", "Diam1", "Diam2", "TankID", "Temperature", "pH", "Treatment") %>%
    filter(TankID != "8t" ) %>% 
      #remove urchin number 8 which died halfway through
    drop_na(Diam1, Diam2) %>% 
    gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>% 
      #gather Diam1 and Diam2 into single column in order to make number of levels of each grouping factor < number of observations (as per error message previously recieved)
    filter(Day == "-13") %>% 
      #sizes on day -24  
    select("TankID", "Diameter", "Temperature", "pH", "Treatment")


AcclimatedMean <- mean(AcclimatedDiameter$Diameter)
AcclimatedSE <- std.error(AcclimatedDiameter$Diameter)
Results of linear mixed model effect on diameter (day -13)
AcclimatedMod <- lmer(Diameter ~ Temperature*pH + (1|TankID), data=AcclimatedDiameter)
anova(AcclimatedMod, type = 2)
## Type II Analysis of Variance Table with Satterthwaite's method
##                  Sum Sq  Mean Sq NumDF DenDF F value  Pr(>F)  
## Temperature    0.015013 0.015013     1    19  0.2154 0.64782  
## pH             0.254996 0.254996     1    19  3.6588 0.07097 .
## Temperature:pH 0.177404 0.177404     1    19  2.5455 0.12711  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions of diameter model (day -13)
#check assumptions
# 1. Normality of residuals
qqPlot(residuals(AcclimatedMod)) #not super normal but... see below

## [1] 42  2
shapiro.test(residuals(AcclimatedMod)) #passes
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(AcclimatedMod)
## W = 0.9577, p-value = 0.09326
hist(residuals(AcclimatedMod)) #looks normalish

# 2. Equal variances
leveneTest(residuals(AcclimatedMod)~InitialDiameter$Treatment) #Passes
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  1.3961 0.2573
##       42
plot(fitted(AcclimatedMod),resid(AcclimatedMod,type="pearson"),col="blue")

Day 1

Assemble Data for diameters on day 1 of experiment when desired conditions were reached (see code)
##After ramp up and desired conditions are reached to begin experiment = Day 1

Day1Diameter <-
  #seperate out the initial sizes of urchins on day 1, when desired conditions were met.
  results %>% 
    select("Day", "Diam1", "Diam2", "TankID", "Temperature", "pH", "Treatment") %>%
    filter(TankID != "8t" ) %>% 
      #remove urchin number 8 which died halfway through
    drop_na(Diam1, Diam2) %>% 
    gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>% 
      #gather Diam1 and Diam2 into single column in order to make number of levels of each grouping factor < number of observations (as per error message previously recieved)
    filter(Day == "1") %>% 
      #sizes on day 1  
    select("TankID", "Diameter", "Temperature", "pH", "Treatment")


Day1Mean <- mean(Day1Diameter$Diameter)
Day1SE <- std.error(Day1Diameter$Diameter)
Results of linear mixed model effect on diameter (day 1)
Day1Mod <- lmer(Diameter ~ Temperature*pH + (1|TankID), data=Day1Diameter)
anova(Day1Mod, type = 2)
## Type II Analysis of Variance Table with Satterthwaite's method
##                 Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## Temperature    0.04177 0.04177     1    19  0.1922 0.66606  
## pH             0.38094 0.38094     1    19  1.7525 0.20128  
## Temperature:pH 0.89969 0.89969     1    19  4.1389 0.05611 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions of diameter model (Day 1)
#check assumptions
# 1. Normality of residuals
qqPlot(residuals(AcclimatedMod)) #not super normal but... see below

## [1] 42  2
shapiro.test(residuals(AcclimatedMod)) #passes
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(AcclimatedMod)
## W = 0.9577, p-value = 0.09326
hist(residuals(AcclimatedMod)) #looks normalish

# 2. Equal variances
leveneTest(residuals(AcclimatedMod)~InitialDiameter$Treatment) #Passes
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  1.3961 0.2573
##       42
plot(fitted(AcclimatedMod),resid(AcclimatedMod,type="pearson"),col="blue")

Growth

Initial as Day -13

Assemble data and calculate means using day -13 as initial diameter (see code)
### GROWTH MODEL 2: If run the same model using day -13 as the first technical day of the experiment. This is the day conditions began to ramp up, so urchins were experiencing gradual increase of stressors.


Initial1 <-
    #seperate out the initial sizes of urchins on day -14, when conditions began to ramp up.
  results %>% 
    select("Day", "Treatment", "Temperature", "pH", "Diam1", "Diam2", "TankID") %>%
    filter(TankID != "8t" ) %>% 
        #remove urchin number 8 which died halfway through
    drop_na(Diam1, Diam2) %>% 
    gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>%
      #gather Diam1 and Diam2 into single column
    filter(Day == "-13") %>%  
      #filter for day at -13 to use this as the intial size.
    select("Temperature", "pH", "Diameter", "TankID")

  mean(Initial1$Diameter)
  std.error(Initial1$Diameter)


Final1 <- 
      #seperate out the final sizes of urchins on day 126
   results %>% 
    select("Day", "Treatment", "Temperature", "pH", "Diam1", "Diam2", "TankID") %>% 
    drop_na(Diam1, Diam2) %>% 
    gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>%
      #gather Diam1 and Diam2 into single column
    filter(Day == "126") %>% 
      #filter for size on last day of experiment
    select("Temperature", "pH", "Diameter","TankID")


Growth1 <- 
  #create column that is growth
  bind_cols(Initial1, Final1) %>% 
    #binds initial (Diameter) and final (Diameter1) sizes to calculate growth
  mutate(Growth = ((Diameter1 - Diameter)/Diameter*100)) %>%
    #add column for growth calculation
  select("Growth")


resultsGrow1 <-
  #table of initial and final size data
  results %>% 
    select("Day", "Treatment", "Temperature", "pH", "TankID", "Diam1", "Diam2") %>% 
    drop_na(Diam1, Diam2) %>% 
    gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>% 
    filter(Day == "126")
resultsGrow1 <- bind_cols(resultsGrow1, Growth1)
  #combines table of initial and final with growth information


#Figure out means of each treatment group
resultsGrow1 %>% 
  dplyr::group_by(Temperature, pH) %>% 
  dplyr::summarise(mean = mean(Growth), s.e. = se(Growth))
Results with linear mixed model effect on Growth (Day -13 as initial).
#LMM for growth:
modGrow1 <- 
  resultsGrow1 %>% 
  lmer(Growth~ Temperature * pH + (1|TankID), data=.)


summary(modGrow1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Growth ~ Temperature * pH + (1 | TankID)
##    Data: .
## 
## REML criterion at convergence: 451.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.97009 -0.36352 -0.06635  0.38074  2.05538 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  TankID   (Intercept) 4120     64.19   
##  Residual              667     25.83   
## Number of obs: 46, groups:  TankID, 23
## 
## Fixed effects:
##                             Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)                  556.515     27.244  19.000  20.427 2.17e-14
## Temperatureheated             56.177     38.529  19.000   1.458    0.161
## pHambient                     -3.451     38.529  19.000  -0.090    0.930
## Temperatureheated:pHambient    9.115     55.833  19.000   0.163    0.872
##                                
## (Intercept)                 ***
## Temperatureheated              
## pHambient                      
## Temperatureheated:pHambient    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Tmprtr pHmbnt
## Tempertrhtd -0.707              
## pHambient   -0.707  0.500       
## Tmprtrhtd:H  0.488 -0.690 -0.690
anova(modGrow1, type="II") #anova{stats}.Compute analysis of variance (or deviance) tables for one or more fitted model objects. Produces Type III Analysis of Variance Table with Satterthwaite's method.
## Type II Analysis of Variance Table with Satterthwaite's method
##                 Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## Temperature    3141.63 3141.63     1    19  4.7100 0.04287 *
## pH                0.68    0.68     1    19  0.0010 0.97489  
## Temperature:pH   17.78   17.78     1    19  0.0267 0.87205  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions of growth model (day -13 as initial)
## ASSSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modGrow1)) 

## [1]  2 27
shapiro.test(residuals(modGrow1)) #passes 
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modGrow1)
## W = 0.95274, p-value = 0.05986
hist(residuals(modGrow1)) #looks normal

# 2. Equal variances
leveneTest(residuals(modGrow1)~resultsGrow1$Treatment) #doesn't pass...
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)  
## group  3  2.9964 0.0413 *
##       42                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(fitted(modGrow1),resid(modGrow1,type="pearson"),col="blue")

Display exploratory interaction plot for growth model (day -13 as initial).
#Interaction plot shows interaction
interaction.plot(resultsGrow1$Temperature,resultsGrow1$pH,resultsGrow1$Growth,
                 type = c("l","p", "b", "o","c"), legend = TRUE,
                 trace.label = deparse(substitute(pH)),
                 fixed = FALSE,
                 xlab = deparse(substitute(Temperature)),
                 ylab = deparse(substitute(Growth)),)

Day 1

Assemble data and calculate means using day 1 as initial diameter and calculate means (see code)
### GROWTH MODEL 1: using day 1 as the first official day of experiment when desired environmental conditions were met (24 days after getting urchins. In that 24 days, they were acclimated and then conditions ramped up)

Initial <-
  #seperate out the initial sizes of urchins on day 1, when desired conditions were met to standardize growth to this measurement.
  results %>% 
    select("Day", "Treatment", "Temperature", "pH", "Diam1", "Diam2", "TankID") %>%
    filter(TankID != "8t" ) %>% 
      #remove urchin number 8 which died halfway through
    drop_na(Diam1, Diam2) %>% 
    gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>% 
      #gather Diam1 and Diam2 into single column in order to make number of levels of each grouping factor < number of observations (as per error message previously recieved)
    filter(Day == "1") %>% 
      #sizes on day 1  
    select("Temperature", "pH", "Diameter", "TankID")

#Mean diameters on day 1
  mean(Initial$Diameter)
## [1] 16.03
  std.error(Initial$Diameter)
## [1] 0.2550688
Final <- 
  #seperate out the final sizes of urchins on day 126
   results %>% 
    select("Day", "Treatment", "Temperature", "pH", "Diam1", "Diam2", "TankID") %>% 
    drop_na(Diam1, Diam2) %>% 
    gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>% 
      #gather Diam1 and Diam2 into single column
    filter(Day == "126") %>% 
      #filter for sizes on day 126
    select("Temperature", "pH", "Diameter","TankID")


Growth <- 
  #create column that is growth
  bind_cols(Initial, Final) %>% 
    #binds initial (Diameter) and final (Diameter1) sizes to calculate growth
  mutate(Growth = ((Diameter1 - Diameter)/Diameter*100)) %>% 
    #add column for growth calculation
  select("Growth")


resultsGrow <-
  #table of initial and final size data
  results %>% 
    select("Day", "Treatment", "Temperature", "pH", "TankID", "Diam1", "Diam2") %>% 
    drop_na(Diam1, Diam2) %>% 
    gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>% 
    filter(Day == "126")
resultsGrow <- bind_cols(resultsGrow, Growth)
#combines table of initial and final with growth information

#Figure out means of each treatment group
resultsGrow %>% 
  dplyr::group_by(Temperature, pH) %>% 
  dplyr::summarise(mean = mean(Growth), s.e. = se(Growth))
## # A tibble: 4 x 4
## # Groups:   Temperature [2]
##   Temperature pH         mean  s.e.
##   <fct>       <fct>     <dbl> <dbl>
## 1 ambient     acidified  327. 16.3 
## 2 ambient     ambient    323. 11.0 
## 3 heated      acidified  352. 15.0 
## 4 heated      ambient    376.  7.15
Growth linear mixed model results (Day 1 as initial)
#LMM for growth:
modGrow <-
  resultsGrow %>% 
  lmer(Growth~ Temperature * pH + (1|TankID), data=.)

#Need to decide on model type that we want to use. 

summary(modGrow)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Growth ~ Temperature * pH + (1 | TankID)
##    Data: .
## 
## REML criterion at convergence: 417.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.90063 -0.32406  0.08141  0.32903  1.94882 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  TankID   (Intercept) 1929.2   43.92   
##  Residual              283.1   16.83   
## Number of obs: 46, groups:  TankID, 23
## 
## Fixed effects:
##                             Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)                  326.595     18.578  19.000  17.580 3.28e-13
## Temperatureheated             25.439     26.273  19.000   0.968    0.345
## pHambient                     -3.424     26.273  19.000  -0.130    0.898
## Temperatureheated:pHambient   27.731     38.073  19.000   0.728    0.475
##                                
## (Intercept)                 ***
## Temperatureheated              
## pHambient                      
## Temperatureheated:pHambient    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Tmprtr pHmbnt
## Tempertrhtd -0.707              
## pHambient   -0.707  0.500       
## Tmprtrhtd:H  0.488 -0.690 -0.690
anova(modGrow, type="II") #anova{stats}.Compute analysis of variance (or deviance) tables for one or more fitted model objects. Produces Type II Analysis of Variance Table with Satterthwaite's method.
## Type II Analysis of Variance Table with Satterthwaite's method
##                 Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## Temperature    1169.41 1169.41     1    19  4.1304 0.05634 .
## pH               74.92   74.92     1    19  0.2646 0.61290  
## Temperature:pH  150.20  150.20     1    19  0.5305 0.47527  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions for growth model (Day 1 as initial)
## ASSSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modGrow)) 

## [1] 43 20
shapiro.test(residuals(modGrow)) #passes
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modGrow)
## W = 0.95491, p-value = 0.07263
hist(residuals(modGrow)) #looks normal

# 2. Equal variances
leveneTest(residuals(modGrow)~resultsGrow$Treatment) #Passes
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  2.1906 0.1033
##       42
plot(fitted(modGrow),resid(modGrow,type="pearson"),col="blue")

Display exploratory interaction plot for growth model (Day 1 as initial)
#Interaction plot shows slight interaction
interaction.plot(resultsGrow$Temperature,resultsGrow$pH,resultsGrow$Growth,
                 type = c("l","p", "b", "o","c"), legend = TRUE,
                 trace.label = deparse(substitute(pH)),
                 fixed = FALSE,
                 xlab = deparse(substitute(Temperature)),
                 ylab = deparse(substitute(Growth)),)

Calcification

Distal End (Tips of spines)

Assemble data for chosen spine tips
### Model 1: calcification at the tips of spines

resultsTip <- 
  #create data using only the SEM images of the tips of spines.
  results %>% 
  select ("Treatment", "Temperature", "pH", "PartOfSpine", "Chosen", "RatioSEM", "TankID") %>% 
    #selects for desired variables. PartOfSpine is tip or base, Chosen indicates which side of the image (arbitrary) is selected for (aka has no dust)
  drop_na(RatioSEM) %>% 
  filter(Chosen == "yes", PartOfSpine == "Tip")
    #filters only for tips and chosen side. No dusty images to be analyzed.
Results for spine tips with linear mixed model
#LMM for calcification at the tips of spines 
modTip <-
  resultsTip %>% 
  lmer(RatioSEM ~ Temperature * pH + (1|TankID), data = .)
## boundary (singular) fit: see ?isSingular
summary(modTip)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: RatioSEM ~ Temperature * pH + (1 | TankID)
##    Data: .
## 
## REML criterion at convergence: 64.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.76413 -0.74507 -0.02052  0.63045  2.07630 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  TankID   (Intercept) 0.0000   0.0000  
##  Residual             0.1356   0.3682  
## Number of obs: 67, groups:  TankID, 23
## 
## Fixed effects:
##                             Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)                  1.69118    0.08930 63.00000  18.938   <2e-16
## Temperatureheated           -0.06662    0.12453 63.00000  -0.535   0.5945
## pHambient                   -0.07462    0.12453 63.00000  -0.599   0.5512
## Temperatureheated:pHambient  0.30557    0.18089 63.00000   1.689   0.0961
##                                
## (Intercept)                 ***
## Temperatureheated              
## pHambient                      
## Temperatureheated:pHambient .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Tmprtr pHmbnt
## Tempertrhtd -0.717              
## pHambient   -0.717  0.514       
## Tmprtrhtd:H  0.494 -0.688 -0.688
## convergence code: 0
## boundary (singular) fit: see ?isSingular
anova(modTip, type="II") #anova{stats}.Compute analysis of variance (or deviance) tables for one or more fitted model objects. Produces Type II Analysis of Variance Table with Satterthwaite's method.
## Type II Analysis of Variance Table with Satterthwaite's method
##                 Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## Temperature    0.10158 0.10158     1    63  0.7492 0.39000  
## pH             0.08185 0.08185     1    63  0.6038 0.44006  
## Temperature:pH 0.38685 0.38685     1    63  2.8534 0.09612 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions for calcification at tip model.

## ASSSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modTip)) #Normal

## [1] 38  8
shapiro.test(residuals(modTip)) #Pass
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modTip)
## W = 0.97737, p-value = 0.2605
hist(residuals(modTip))

# 2. Equal variances
leveneTest(residuals(modTip)~resultsTip$Treatment) #Pass
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  0.4545  0.715
##       63
plot(fitted(modTip),resid(modTip,type="pearson"),col="blue")

Conduct post hoc analysis on calcification at tip.

## POST HOC ANALYSIS 
emmTip <- emmeans(modTip, ~Temperature*pH, adjust = "tukey")  #Estimated marginal means (Least-squares means)
multcomp::cld(emmTip) #create a compact letter display of all pair-wise comparisons
##  Temperature pH        emmean     SE   df lower.CL upper.CL .group
##  ambient     ambient     1.62 0.0868 18.2     1.43     1.80  1    
##  heated      acidified   1.62 0.0868 18.2     1.44     1.81  1    
##  ambient     acidified   1.69 0.0900 18.2     1.50     1.88  1    
##  heated      ambient     1.86 0.0993 16.6     1.65     2.07  1    
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
pwpp(emmTip) #Constructs a plot of P values associated with pairwise comparisons of estimated marginal means.

Display interaction plot for exploration of calcification at tip.

#Interaction plot for cacification at the tip of the spines - interaction
interaction.plot(resultsTip$Temperature,resultsTip$pH,resultsTip$RatioSEM,
                 type = c("l","p", "b", "o","c"), legend = TRUE,
                 trace.label = deparse(substitute(pH)),
                 fixed = FALSE,
                 xlab = deparse(substitute(Temperature)),
                 ylab = deparse(substitute(Growth)),)

Proximal End (Base of spines)

Assemble data for chosen spine bases

### Model 2: calcification in the base of the spines

resultsBase <-
   #create data using only the SEM images of the base of spines.
  results %>% 
  select ("Treatment", "Temperature", "pH", "PartOfSpine", "Chosen", "RatioSEM", "TankID") %>% 
    #selects for desired variables. PartOfSpine is tip or base, Chosen indicates which side of the image (arbitrary) is selected for (aka has no dust)
  drop_na(RatioSEM) %>% 
  filter(Chosen == "yes", PartOfSpine == "Base")
    #filters only for tips and chosen side. No dusty images to be analyzed.

Results for spine bases with linear mixed model

#LMM for calcification at the tips of spines 
modBase <- 
  resultsBase %>% 
  lmer(RatioSEM~ Temperature * pH + (1|TankID), data=.)

summary(modBase) 
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: RatioSEM ~ Temperature * pH + (1 | TankID)
##    Data: .
## 
## REML criterion at convergence: 98.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5535 -0.6315 -0.1529  0.6918  1.9937 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  TankID   (Intercept) 0.1303   0.3610  
##  Residual             0.1597   0.3996  
## Number of obs: 68, groups:  TankID, 23
## 
## Fixed effects:
##                             Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)                   2.1440     0.1749 18.8077  12.257 2.06e-10
## Temperatureheated             0.3566     0.2487 19.1707   1.434   0.1677
## pHambient                     0.6524     0.2474 18.8077   2.637   0.0163
## Temperatureheated:pHambient  -0.2224     0.3594 18.9804  -0.619   0.5434
##                                
## (Intercept)                 ***
## Temperatureheated              
## pHambient                   *  
## Temperatureheated:pHambient    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Tmprtr pHmbnt
## Tempertrhtd -0.703              
## pHambient   -0.707  0.497       
## Tmprtrhtd:H  0.487 -0.692 -0.688
anova(modBase, type="II") #anova{stats}.Compute analysis of variance (or deviance) tables for one or more fitted model objects. Produces Type III Analysis of Variance Table with Satterthwaite's method.
## Type II Analysis of Variance Table with Satterthwaite's method
##                 Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## Temperature    0.30743 0.30743     1 18.992  1.9251 0.181362   
## pH             1.48873 1.48873     1 18.960  9.3223 0.006553 **
## Temperature:pH 0.06114 0.06114     1 18.980  0.3829 0.543424   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions for calcification at base model.

## ASSSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modBase)) #Normal

## [1] 33 37
shapiro.test(residuals(modBase)) # Pass
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modBase)
## W = 0.97304, p-value = 0.1472
hist(residuals(modBase))

# 2. Equal variances
leveneTest(residuals(modBase)~resultsBase$Treatment) #Pass
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3   1.068  0.369
##       64
plot(fitted(modBase),resid(modBase,type="pearson"),col="blue")

Conduct post hoc analysis on calcification at base

## POST HOC ANALYSIS 
emmBase <- emmeans(modBase, ~Temperature*pH, adjust = "tukey")
  #Estimated marginal means (Least-squares means)
multcomp::cld(emmBase)
##  Temperature pH        emmean    SE   df lower.CL upper.CL .group
##  ambient     acidified   2.14 0.175 18.8     1.78     2.51  1    
##  heated      acidified   2.50 0.177 19.5     2.13     2.87  12   
##  ambient     ambient     2.80 0.175 18.8     2.43     3.16  12   
##  heated      ambient     2.93 0.192 18.8     2.53     3.33   2   
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
  #create a compact letter display of all pair-wise comparisons
pwpp(emmBase)

  #Constructs a plot of P values associated with pairwise comparisons of estimated marginal means.

Display interaction plot for exploration of calcification at base

#Interaction plot for calcification at the base of the spines - no interaction
interaction.plot(resultsBase$Temperature,resultsBase$pH,resultsBase$RatioSEM,
                 type = c("l","p", "b", "o","c"), legend = TRUE,
                 trace.label = deparse(substitute(pH)),
                 fixed = FALSE,
                 xlab = deparse(substitute(Temperature)),
                 ylab = deparse(substitute(RatioSEM)),)

Dropped Spines

Assemble data for dropped spines
## Number of loose spines counted at the bottom of tanks at the end of the experiment. These were either broken or shed. Note: these spines were not collected, just counted, so can not determine whether they were broken or shed.

resultsDropped <-
  #create data using only the needed variables.
  results %>% 
  select ("Treatment", "Temperature", "pH", "SpineCount", "TankID") %>% 
  drop_na(SpineCount)
Analyse dropped spines with linear mixed effect model.
##LM of dropped spines - can't make it LMM because error: number of levels of each grouping factor must be < number of observations - only have one count for each TankID.
modDropped <- 
  resultsDropped %>% 
  lm(SpineCount~ Temperature * pH, data=.)

summary(modDropped)
## 
## Call:
## lm(formula = SpineCount ~ Temperature * pH, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.0000  -6.0833  -0.1667   0.4000  20.0000 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   21.167      3.640   5.814 1.33e-05 ***
## Temperatureheated             -1.167      5.148  -0.227  0.82315    
## pHambient                    -21.000      5.148  -4.079  0.00064 ***
## Temperatureheated:pHambient    1.600      7.461   0.214  0.83248    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.917 on 19 degrees of freedom
## Multiple R-squared:  0.6088, Adjusted R-squared:  0.547 
## F-statistic: 9.855 on 3 and 19 DF,  p-value: 0.0003899
anova(modDropped) #anova{stats}.Compute analysis of variance (or deviance) tables for one or more fitted model objects. Produces Type III Analysis of Variance Table with Satterthwaite's method.
## Analysis of Variance Table
## 
## Response: SpineCount
##                Df  Sum Sq Mean Sq F value    Pr(>F)    
## Temperature     1    1.52    1.52  0.0192    0.8914    
## pH              1 2345.78 2345.78 29.4995 3.062e-05 ***
## Temperature:pH  1    3.66    3.66  0.0460    0.8325    
## Residuals      19 1510.87   79.52                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions for model for dropped spines
##ASSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modDropped)) #not normal

## 1458 1460 
##    2    4
shapiro.test(residuals(modDropped)) #faaaiiilll
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modDropped)
## W = 0.86706, p-value = 0.005636
hist(residuals(modDropped))

# 2. Equal variances
leveneTest(residuals(modDropped)~resultsDropped$Treatment) # barely pass
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  3  2.8487 0.06483 .
##       19                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(fitted(modDropped),resid(modDropped,type="pearson"),col="blue")

Analysis for dropped spines does not meet assumptions. Attempt a square root transformation.
##DOES NOT MEET ASSUMPTIONS so ATTEMPT AT Transformation:
resultsDropped$tdata<-(resultsDropped$SpineCount)^1/2
  # transformation
modDropped1 <- lm(tdata~ Temperature * pH, data=resultsDropped)
  # run lm model of transformed data

###ASSUMPTIONS for Transformation
# 1. Normality of residuals
qqPlot(residuals(modDropped1)) #not normal

## 1458 1460 
##    2    4
shapiro.test(residuals(modDropped1)) #fail
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modDropped1)
## W = 0.86706, p-value = 0.005636
hist(residuals(modDropped1))

# 2. equal variances
bartlett.test(residuals(modDropped1)~resultsDropped$Treatment) #fail
## 
##  Bartlett test of homogeneity of variances
## 
## data:  residuals(modDropped1) by resultsDropped$Treatment
## Bartlett's K-squared = 43.282, df = 3, p-value = 2.144e-09
plot(fitted(modDropped1),resid(modDropped1,type="pearson"),col="blue")

summary(modDropped1)
## 
## Call:
## lm(formula = tdata ~ Temperature * pH, data = resultsDropped)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.0000 -3.0417 -0.0833  0.2000 10.0000 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  10.5833     1.8202   5.814 1.33e-05 ***
## Temperatureheated            -0.5833     2.5742  -0.227  0.82315    
## pHambient                   -10.5000     2.5742  -4.079  0.00064 ***
## Temperatureheated:pHambient   0.8000     3.7304   0.214  0.83248    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.459 on 19 degrees of freedom
## Multiple R-squared:  0.6088, Adjusted R-squared:  0.547 
## F-statistic: 9.855 on 3 and 19 DF,  p-value: 0.0003899
Anova(modDropped1, type=2) #test for significance of model
## Anova Table (Type II tests)
## 
## Response: tdata
##                Sum Sq Df F value    Pr(>F)    
## Temperature      0.23  1  0.0118    0.9146    
## pH             586.44  1 29.4995 3.062e-05 ***
## Temperature:pH   0.91  1  0.0460    0.8325    
## Residuals      377.72 19                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Transfomraiton (of any usual kind) of dropped spines was not successful, so we use a non parametric Kruskal Wallis Test.
### Transformations not work, so use non parametric:
kruskal.test(resultsDropped$SpineCount~resultsDropped$Treatment)# SIG*** (p=0.0005563)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  resultsDropped$SpineCount by resultsDropped$Treatment
## Kruskal-Wallis chi-squared = 17.505, df = 3, p-value = 0.0005563
Conduct Dunn Post Hoc Test for dropped spines
## POST HOC Pairwise
dunn <- 
  resultsDropped %>% 
  dunnTest(SpineCount ~ Treatment,
           method = "holm", kw = TRUE, data = .)
dunnph <- dunn$res
cldList(P.adj ~ Comparison, data = dunnph, threshold = 0.05)
##   Group Letter MonoLetter
## 1   AMB      a        a  
## 2    OA      b         b 
## 3     T     ac        a c
## 4  T/OA     bc         bc
Display interaction plot for dropped spines
#Interaction plot for dropped spines - no interaction
interaction.plot(resultsDropped$Temperature,resultsDropped$pH,resultsDropped$SpineCount,
                 type = c("l","p", "b", "o","c"), legend = TRUE,
                 trace.label = deparse(substitute(pH)),
                 fixed = FALSE,
                 xlab = deparse(substitute(Temperature)),
                 ylab = deparse(substitute(SpineCount)),)